SU(2) Landau gluon propagator on a 140 3 lattice 
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We present a numerical study of the gluon propagator in 
lattice Landau gauge for three-dimensional pme-SU(2) lat- 
tice gauge theory at couplings /3 = 4.2, 5.0, 6.0 and for lattice 
volumes V = 40 3 , 80 3 , 140 3 . In the limit of large V we ob- 
serve a decreasing gluon propagator for momenta smaller than 
Pdec = 350t5°° MeV. Data are well fitted by Gribov-like for- 
mulae and seem to indicate an infra-red critical exponent k, 
slightly above 0.6, in agreement with recent analytic results. 
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I. INTRODUCTION 

The study of the infra-red (IR) limit of QCD is of 
central importance for understanding the mechanism of 
quark confinement and the dynamics of partons at low 
energy. Despite being non-gauge- invariant, the gluon 
propagator is a powerful tool in this (non-perturbative) 
investigation [1]. In particular, it would be interesting to 
express it in a closed form for recovering the phenomenol- 
ogy of Pomeron exchange from first principles [2] . 

Studies of the coupled set of Dyson-Schwinger equa- 
tions for gluon and ghost propagators in Landau gauge 
predict for the gluon propagator an IR behavior of the 
form D(p) — p 4K - 2 [implying D(0) — if K > 0.5]. 
The available predictions for the IR exponent are n € 
[0.52, 1.00] in the four-dimensional case [3,4] and k w 
0.648 or k — 0.75 in three dimensions [4]. 

Furthermore, in the minimal Landau gauge, the gauge- 
fixed configurations belong to the region f2 of transverse 
configurations, for which the Faddeev-Popov operator is 
non-negative. This implies a rigorous inequality [5] for 
the Fourier components of the gluon field A^(x) and a 
strong suppression of the gluon propagator in the IR 
limit. In particular, for dimension d and infinite volume, 
it is proven that the (unrenormalized) gluon propagator 
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is less singular than p 2 ~ d and that, very likely, it van- 
ishes in the IR limit [5]. A vanishing gluon propagator 
at p = 0, given by the form p 2 /(p 4 + A 4 ), was also ob- 
tained by Gribov [6]. Here the mass scale A arises when 
the configuration space is restricted to the region Q. A 
generalization of this expression has been introduced in 
Ref. [7] as an Ansatz for a non-perturbative solution of 
the gluon Dyson-Schwinger equation. 

Numerical studies [8,9] have now established that 
the gluon propagator in lattice Landau gauge shows a 
turnover in the IR region and attains a finite value for 
p = 0. Evidences of a decreasing gluon propagator for 
small p have been obtained in the 4d SU(2) and SU(3) 
cases (but only in the strong-coupling regime) [10,11] in 
the 3d SU(2) case (also in the scaling region) [8,12,13], 
in the 3d SU(2) adjoint Higgs model [13], in the 4c? 
SU(2) case at finite temperature [14] and for the equal- 
time three-dimensional transverse gluon propagator in 4d 
SU{2) Coulomb gauge [15]. In this last case, one obtains 
an excellent fit of the transverse propagator by a Gribov- 
like formula. 

This work aims to verify the possibility of using 
Gribov-like formulae to fit data of the gluon propaga- 
tor also in Landau gauge. At the same time, we will 
try to obtain a value for the IR critical exponent k to 
be compared to the analytic determinations mentioned 
above. In order to probe the infinite-volume limit and 
the IR region we consider the three-dimensional case and 
the SU (2) group, using lattice sizes up to 140 3 . Note that 
the study of the gluon propagator in three dimensions is 
also of interest in finite-temperature QCD [16]. 

II. NUMERICAL SIMULATIONS 

We consider the standard Wilson action for SU(2) 
lattice gauge theory in three dimensions with periodic 
boundary conditions. The numerical code is entirely par- 
allelized using MPI. (Technical details and a study of the 
code performance are left for a subsequent work [17].) 
For the construction of staples we follow Ref. [18]. For 
the random number generator we use a double-precision 
implementation of RANLUX (version 2.1) with luxury 
level set to 2. Computations were performed on the PC 



cluster at the IFSC-USP. The system has 16 nodes and 
a server with 866 MHz Pentium III CPU and 256/512 
MB RAM memory. The machines are connected with a 
100 Mbps full-duplex network. The total computer time 
used for the runs was about 80 days on the full cluster. 

In Table I we report, for each coupling (3 and lattice 
volume V , the parameters used for the simulations. All 
our runs start with a random gauge configuration. For 
thermalization we use a hybrid overrelaxed (HOR) algo- 
rithm [19]. Each HOR iteration consists of one heat- 
bath sweep over the lattice followed by to micro-canonical 
sweeps. We did not try to find the best tuning for to; we 
use to = 4 for V = 40 3 , 80 3 and m = 5 for V = 140 3 . 
In order to optimize the heat-bath code, we implement 
two different SU{2) generators, namely methods 1 and 2 
described in [20, Appendix A] with h cuto ff = 2. 

For the numerical gauge fixing we use the stochastic 
overrelaxation (SOR) algorithm [21,22] with even/odd 
update. In Table I we report the value of the tuning 
parameter p sor used for each pair {(3,V). We stop the 
gauge fixing when the average value of [(V • A) b (x)] 2 is 
smaller than 10~ 12 . (For a definition of the lattice gauge 
field A^x) and of the lattice divergence V we refer to [8].) 
We note that half of the configurations for V = 80 3 were 
done using the so-called Cornell method [21,22], with tun- 
ing parameters a corn = 0.325, 0.32, 0.316 respectively at 
/3 = 4.2,5.0,6.0. In fact, the Cornell method is somewhat 
faster than the SOR algorithm and leads to a gauge fix- 
ing of comparable quality if one uses an even/odd update 
[22] . A good estimator of the quality of the gauge fixing is 
the quantity £q (see eq. 6.8 in Ref. [22]), which should be 
zero when the configuration is gauge-fixed. By averaging 
over the gauge-fixed configurations, we find (at (3 = 4.2) 
that the ratio between the final and the initial values of 
Eq is (in 95% of the cases) about 5.3 x 10~ 10 with the 
Cornell method and 2.4 x 10~ n with the SOR algorithm. 
At the same time, the average CPU-time needed for up- 
dating each site variable is about 11% smaller for the 
Cornell method. In any case, the CPU-time for gauge 
fixing was quite significant for the large lattices. In order 
to go to even larger lattices one should probably imple- 
ment a global algorithm such as the Fourier acceleration 
method [21,22] with the multigrid or conjugate gradient 
implementations introduced in Ref. [23] , which are highly 
parallelizable. 

We ran the 40 3 lattices on a single node, the 80 3 lat- 
tices on two nodes and the 140 3 lattices on four nodes. 
The parallclization of the code worked well. In fact, for 
runs on 2 and 4 nodes, we obtain speed-up factors 1.82 
and 3.41 for the heat-bath link update. For the micro- 
canonical link update the factors are respectively 1.87 
and 3.77 and for the SOR site update we get 1.90 and 
3.72. 

For each (3 we evaluate the average plaquette (1^1,1) 
(see Table II). Results are in agreement with the data 
reported in Ref. [8], but we now have smaller statisti- 
cal errors. (The data have been analyzed using various 



methods, described in footnote 4 of Ref. [8]. Here we 
always report the largest error found.) We also evalu- 
ate the tadpole-improved coupling (3i = f3(Wi t i). In 
this way, by using the fit given in eq. 2 and Table IV 
of Ref. [24] we calculate the string tension ^fo in lattice 
units (see Table II) and the inverse lattice spacing a^ 1 
using the input value ^/a = 0.44 GeV. The fit is valid 
for (3 > 3.0, i.e. the couplings (3 considered here are well 
above the strong-coupling region. Let us notice that, if 
we compare the data for the string tension (in lattice 
units) with data obtained for the SU (2) group in four 
dimensions [25, Table 3], then our values of (3 correspond 
to (3 ~ 2.28,2.345,2.41 in the four dimensional case. Fi- 
nally, in the same table we report the lattice spacing a in 
fm and the smallest non-zero momentum (in McV) that 
can be considered for each (3. Thus, with the lattice vol- 
umes and the (3 values used here we are able to consider 
momenta as small as 59 MeV (in the deep IR region) and 
physical lattice sides almost as large as 25 fm. 

In this work we did not do a systematic study of 
Gribov-copy effects for the gluon propagator. However, 
we compared data obtained using the SOR and the Cor- 
nell gauge-fixing methods. In principle, different gauge- 
fixing algorithms — or even the same algorithm with dif- 
ferent values of the tuning parameter [26] — can generate 
different Gribov copies starting from the same thermal- 
ized configuration. Thus, this comparison provides an 
estimate of the bias (Gribov noise) introduced by the 
gauge-fixing procedure. We found that, in most cases, 
the difference between the two sets of gluon-propagator 
data is within 1 standard deviation and that, in all cases, 
it is smaller than 2 standard deviations. Moreover, this 
difference did not show any systematic effect, suggesting 
that the influence of Gribov copies on the gluon prop- 
agator (if present) is of the order of magnitude of the 
numerical accuracy. This is in agreement with previous 
studies in Landau gauge for the SU (2) and SU(3) groups 
in four dimensions [10,26,27]. A similar result has also 
been obtained for the Coulomb gauge [28]. Note that, in 
the U(l) case [29], Gribov copies can affect the behavior 
of the photon propagator, making it difficult to repro- 
duce the known perturbative behavior in the Coulomb 
phase. The situation is very different for the SU(2) and 
SU(3) cases, at least when considering the lattice Lan- 
dau gauge. In fact, as said in the Introduction, in the 
non-Abelian case the minimal-Landau-gauge condition 
implies [5] the positiveness of the Faddcev-Popov ma- 
trix and a strict bound for the Fourier components of the 
gluon field A^k). The bound applies to all Gribov copies 
obtained with the numerical gauge fixing. Thus, if the 
behavior of the gluon propagator is determined by this 
bound [5,6], then this behavior should be the same for 
all lattice Gribov copies. This would explain why we do 
not see Gribov-copy effects here. Clearly, the same result 
does not apply to the U(l) theory, since in this case the 
Faddeev-Popov matrix is independent of the gauge field. 
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III. RESULTS AND CONCLUSIONS 

We evaluate the lattice gluon propagator D(k) and 
study it as a function of the lattice momentum p 2 (k) 
(see Ref. [8] for definitions). In our simulations we 
consider, for each gauge-fixed configuration, all vectors 
k = (k x ,k y ,kt) with only one component different from 
zero and average over the three directions. For the gluon 
propagator we analyze the data by estimating the statis- 
tical error with three different methods: standard devi- 
ation, jack-knife with single-data elimination and boot- 
strap (with 10000 samples). We found that the results 
obtained are in agreement in all cases. Here we always 
use the standard-deviation error. 

Let us recall that in the three-dimensional case the 
coupling g 2 has dimensions of mass. Thus, in order to 
obtain a dimcnsionless lattice coupling we have to set 
(3 — A/(ag 2 ). Then, with our notation [8], the quan- 
tity aD(k) approaches g 2 D^ cont \k) / A in the continuum 
limit, where D^ cont \k) is the unrenormalized continuum 
gluon propagator. 

In order to compare lattice data at different (3's, we ap- 
ply the matching technique described in [30, Sec. V.B.2]. 
(Note that we have already determined the lattice spac- 
ing a, as described above.) We start by checking for 
finite-size effects, comparing data at different lattice sizes 
and same (3 value. In this way, we find (for each (3) a 
range of ultra-violet (UV) momenta for which the data 
are free from finite- volume corrections. We then perform 
the matching using data for these momenta and V = 40 3 , 
since for this lattice volume the errors are smallest (about 
1%). In particular, when matching data obtained at two 
different values of (3, we first interpolate the data for the 
larger (3 (the fine lattice) using a spline. Finally, we find 
the multiplicative factor Rz = Z(af) / Z(a c ) correspond- 
ing to the best fit of the (multiplied) coarse-lattice data to 
the same spline. The error of Rz is estimated using a pro- 
cedure similar to the one described in [30]. The method 
works very well (see Figure 1). Notice that we did not 
fix the remaining global factor Z imposing a renormal- 
ization condition, as done for example in Ref. [31]. Our 
case is equivalent to setting Z(a) = 1 at f3 = 6.0. 

The data obtained after the matching are shown for 
V = 80 3 in Figure 1. Clearly, we find that the gluon 
propagator decreases in the IR limit for momenta smaller 
than pdec, which corresponds to the mass scale A in a 
Gribov-like propagator. From the plot we can estimate 
p dec = 350±5o° MeV, in agreement with Ref. [8]. 

In Figure 2 we plot the rescaled gluon propagator at 
zero momentum, namely aD(0)/Z(a), as a function of 
the inverse lattice side L^ 1 = l/(aN) in physical units 
(fm _1 ). We see that aD(0)/Z(a) decreases monotoni- 
cally as L increases, in agreement with Ref. [9]. It is in- 
teresting to notice that these data can be well fitted using 
the simple Ansatz d + b/L c both with d = and d ^ 
(see Figure 2). In order to decide for one or the other 
result one should go to significantly larger lattice sizes. 



We plan [32] to extend these simulations to [3 — 3.4 and 
lattice sizes up to 260 3 , allowing us to consider a value 
L _1 w 0.017 frn -1 . (This requires running in parallel on 
all nodes of our PC cluster.) 

Following Ref. [15] we fit the data using a Gribov-like 
(or Stingl-likc) formula 



D(p) = 



zp 



2a 



y 2 + (p 2 + x) 



2 ' 



(1) 



where z, s, a, x and y are fitting parameters. For non- 
negative a this implies a finite gluon propagator in the 
IR limit, with a behavior given by D(p) oc (s + zp 2a ). 
If y 2 > this form corresponds to a propagator with 
poles at m\ = —x ± iy, while if y 2 < the poles are 
real. Finally, note that the IR exponent n considered in 
studies of Dyson-Schwingcr equations is given in terms of 
a by k = (1 + a)/2, assuming D(0) = 0. Results of our 
fits using the un-rescaled lattice data D(k) are reported 
in Table III. We see that a decreases when the physical 
lattice volume increases. Also, we always get y 2 > and 
x < \y\, which seems to support the scenario of purely 
imaginary poles found in [15]. Let us recall that the poles 
of the gluon propagator are gauge-invariant (at all orders 
in perturbation theory) [33]. 

Notice that this fitting form leads in general to the 
wrong UV behavior, name ly D(p) ~ p 2 °--^. This is not 
a serious problem, since the largest momentum we can 
consider is about 3.5 GeV, i.e. we are not really exploring 
the UV limit. On the other hand, the exponent a in 
(1) plays a role both in the IR and in the UV regimes. 
Thus, the values obtained for a probably correspond to 
averaging the behavior of the gluon propagator in these 
two regions. In particular, since we expect D(p) ~ p~ 2 
in the UV limit, it is likely that the IR behavior of the 
propagator be given by a smaller exponent a than the 
ones reported in Table III. To check this we set x = 
and introduce an anomalous dimension 7: 



D(p) 



s + zp- 



2a 



y 2 + p2(i+ 7 ) ' 



(2) 



Results for this fitting form are also reported in Table 
III. Indeed, we get values of a smaller than in the pre- 
vious case (and still decreasing with increasing physical 
lattice volume) . For the anomalous dimension we obtain 
7 Ri 0.65, with small volume- and /3-dcpcndcnce. The 
problem with eq. (2) is that the introduction of 7 com- 
promises the pole interpretation. 
Finally, one can also try the form 



D(p) 



2a 



s + zp- 

(y 2 +P 4 r ' 



(3) 



This corresponds to a propagator with purely imaginary 
poles m\ = Hy and at the same time allows the data 
to select the IR and the UV behaviors separately. The 
problem in this case is that the fit is unstable for small 
lattice volumes. On the contrary, for V = 140 3 the fit 
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works well and we obtain a = 0.27(6), 0.29(7), 0.38(8) 
respectively for (3 = 4.2, 5.0, 6.0 and 7 w 0.72. Hence, 
our a values are of the order of 0.3, corresponding to 
k w 0.65. (Again, there is a decrease of a when the 
physical lattice volume increases.) 

In order to check for possible effects from the breaking 
of rotational invariance [34] we redid our fits substituting 
p 2 (k) by p 2 (k) =p 2 (k) +pW(jfe)/12 (see Ref. [35]). One 
expects this modification to play an important role in the 
UV limit and to have a small effect in the IR case. For 
all (3 values and lattice volumes we obtain good fits and 
results similar to those reported above. In particular, 
using eq. (1) one still finds y 2 > and x < \y\, support- 
ing the scenario of purely imaginary poles. At the same 
time, for the three fitting functions, a decreases when the 
physical lattice volume increases, but its value is about 
20 — 30% larger than the one reported in Table (III) and 
above. In particular, using eq. (3) and data for V = 140 3 
we obtain a = 0.32(7), 0.39(8), 0.59(12) respectively for 
(3 = 4.2, 5.0,6.0 and 7 w 0.7. This implies slightly larger 
values for k = (l + a)/2. The analysis of these discretiza- 
tion effects and their influence on the extrapolation of a 
and k to the continuum limit will be done elsewhere [32]. 

We have confirmed, by numerical simulations in the 
scaling region, that the transverse gluon propagator in 
3d SU(2) Landau gauge is a decreasing function for mo- 
menta p < 350 MeV (attaining a finite value at p = 0) . 
Also, the data are well fitted by the Gribov-likc formulae 
(l)-(3) and we obtain an IR critical exponent k in agree- 
ment with recent analytic results. In order to eliminate 
remaining discretization and finite- volume effects [and in 
particular to check if D(0) = 0], we need to simulate at 
larger lattice volumes and at other values of (3. 
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p/a (GeV) p/a (GeV) 

FIG. 1. Plot of the rescaled gluon propagator as a function of the lattice momentum for V = 80 and f3 — 4.2 (x), 5.0 (□), 
6.0 (O). The second plot shows only the IR region. Error bars are obtained from propagation of errors. 
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FIG. 2. Plot of the rescaled gluon propagator at zero 
momentum as a function of the inverse lattice side for 
(3 = 4.2 (x), 5.0(D), 6.0(O). We also show the fit of the 
data using the Ansatz d + b/L c both with d — and d =fc 0. 
Error bars are obtained from propagation of errors. 

TABLE II. For each coupling j3 we report the value of 
the average plaquette the string tension ^fa in lat- 

tice units, the lattice spacing in fm and the smallest non-zero 
momentum (in MeV) for the lattice volume V = 140 . Error 
bars for have been obtained taking into account the 

integrated autocorrelation time of the HOR algorithm. All 
the other error bars come from propagation of errors. 

P (Wi.i) a (fm) p min (MeV) 

4.2 0.741861(2) 0.387(3) 0.174(1) 59.0(4) 

5.0 0.786869(2) 0.314(2) 0.1407(8) 62.9(4) 

6.0 0.824780(1) 0.254(1) 0.1138(5) 77.8(4) 



TABLE I. The pairs (f3,V) considered for the simulations, 
the number of configurations, the numbers of HOR sweeps 
used for thermalization and between two consecutive config- 
urations (used for evaluating the gluon propagator) and the 
parameter p sor used by the SOR algorithm. 



13 V Configurations Thermalization Sweeps p s , 



4.2 
4.2 
4.2 


40 3 
80 3 
140 3 


400 
200 
30 




1100 
2200 
2750 


100 0.70 
200 0.80 
250 0.88 


5.0 
5.0 
5.0 


40 3 
80 3 
140 3 


400 
200 
30 




1320 
2420 
3080 


120 0.69 
220 0.80 
280 0.88 


6.0 
6.0 
6.0 


40 3 
80 3 
140 3 


400 
200 
30 




1540 
2680 
3300 


140 0.68 
240 0.80 
300 0.87 


TABLE III. Fit of the data using eqs. (1) and (2). In all 
cases x 2 /d.o.f. was of order 1. Note that for a lattice volume 
V = iV 3 we have 1 + N/2 data points and that all points have 
been used for the fits. 


a 


V 




Fit 1 




Fit 2 






a 


X 


\y\ 


a y* 


4.2 
4.2 
4.2 


40 3 
80 3 
140 3 


0.69(2) 
0.66(2) 
0.61(4) 


0.17(2) 
0.17(2) 
0.19(4) 


0.42(1) 
0.36(1) 
0.35(3) 


0.48(4) 0.29(1) 
0.46(4) 0.24(1) 
0.38(6) 0.25(2) 


5.0 
5.0 
5.0 


40 3 
80 3 
140 3 


0.77(3) 
0.71(2) 
0.71(3) 


0.11(2) 
0.11(1) 
0.10(2) 


0.28(2) 
0.214(9) 
0.22(2) 


0.53(6) 0.152(8) 
0.45(3) 0.118(6) 
0.44(6) 0.12(1) 


6.0 
6.0 
6.0 


40 3 
80 3 
140 3 


0.86(3) 
0.84(2) 
0.80(2) 


0.069(8) 
0.032(6) 
0.037(8) 


0.20(2) 
0.166(8) 
0.123(7) 


0.64(6) 0.082(6) 
0.65(5) 0.051(3) 
0.55(6) 0.040(4) 
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